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ABSTRACT 

A method is given for the study of families of periodic motion using 
trigonometric series to represent the individual solutions. The continu- 
ous deformations of a periodic orbit along a family are represented by 
the variations of the trigonometric coefficients with respect to a parameter. 

For the applications the series are truncated, while the coefficients 
and their variations are determined numerically. In this form, the con- 
tinuation with respect to the parameter is given by a mapping f: R n - R n 
of the space of coefficients into itself. The values of the coefficients are 
then improved by another mapping which is a contraction operator in some 
neighborhood of the fixed point representing the solution. 

The method is applied to the natural families defined by Wintner and 
to the families of the first and second kinds of Poincare' in the restricted 
problem of three bodies. 
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PERIODIC ORBITS IN TRIGONOMETRIC SERIES 
INTRODUCTION 

Trigonometric series have been used very extensively in the study of periodic 
orbits and families of periodic motion in the restricted problem of three bodies. 
Many of these studies depend on the presence of small parameters in the problem, 
and the results consist of trigonometric terms whose coefficients are power 
series in the small parameters. That is to say, the results are obtained in the 
form of Poisson series. The same techniques can be applied using numerical 
values for the parameters, and this is referred to as a semi-analytic development. 
However, many families of periodic orbits have not yielded to these methods and 
have thus far been explored only by numerical integration. 

The purpose of the present study is to obtain semi-analytic results for some 
of these more difficult families of periodic orbits and to place the continuation 
of these families on a sound basis when trigonometric series with numerical 
coefficients are used in a method of successive approximations. 

Cases of near and exact resonance are of special interest and are handled 
without difficulty using the present technique. The periodic matrix of the varia- 
tional equations plays an essential role in the method, and the linear stability 
analysis is done in the usual way. 

A good discussion of the method, as applied to isolated periodic solutions 
has been given by Urabe (1965). Numerical examples are given by Urabe and 
Reiter (1966). More recent results and references are given by Stokes (1969). 

The successive approximations are similar to those of Bennett and Palmore 
(1968) in that each step yields a periodic function which is an approximate 
solution of the original equations of motion. The variations of the series 
coefficients with re ^ect to a local parameter give a representation of the 
functions x^t) and y^t) discussed by Depr it and Henrard (1967) (page 160). 

The method described in this study has been applied in the computation of 
periodic orbits in the restricted problem of four bodies by Kolenkiewicz and 
Carpenter (1967 and 1968) and in the restricted problem of three bodies by 
Carpenter and Stumpff (1968). Further applications are being made by Deprit 
and Carpenter for locating elements of many new families or orbits in the 
problem of three bodies. 
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THE ITERATION PROCESS FOR FIXED PERIOD 


The method is applicable to a wide range of problems, and a general formu- 
lation is possible. However, each class of problems has its own interesting 
features, and it is usually possible to greatly improve the efficiency by taking 
advantage of the particular forms of the equations and the solutions. 

The formulation will bo gl.on for the restricted problem of three-bodies 
considering symmetric orbits which lie in the plane of motion of the primaries. 
Furthermore, the method will be developed as a modification of an iterative 
general perturbations technique, so a two-body reference orbit is used as the 
starting point. For convenience the reference orbit is circular with respect to 
the central primary whose mass is put equal to M (referred to as the Sun). 

Let a and n be the radius and mean motion respectively of the circular 
reference orbit of the infinitesimal particle (referred to as the minor planet), 
and let k be the Gaussian constant so that 


n 2 a 3 = k 2 M . 

Let the true position of the minor planet be given by 

r = (1 + a) r Q + /3w 

where r 0 is the position in the reference orbit, 

w _ 1 dr o 
n dt 


is the vector of length a in the direction of the velocity in the reference orbit, 
and the quantities a and /3 will represent the periodic deviations from the 
circular motion. 

Let a' and n' be the radius and mean motion respectively of the circular 
motion of Jupiter, whose mass is demoted by m, so that 

n' 2 a' 3 = k 2 (M + m) or n' 2 a' 3 = k 2 M (1 + m' ) with m' * m/M 


and 


a '/a = (1 + m' ) 1//3 v 2//3 


2 


where 


v = n/n' . 


For the cases where n > n', the relative reference motions will be periodic 
with period 


T 


2 77 


n - n 


I 


and the orbit to be determined will have the period 


with 


T 


2tt 

IT 


<1 

for some positive integer value of The integer l is the winding number or 
index of the orbit with respect to the Sun in the synodic coordinate system. Thus 
the trigonometric argument is 


e = Nt. 


For symmetric orbits the epoch of time is chosen such that the minor planet 
crosses the Sun-Jupiter line at t = 0. 

With respect to the mean anomaly, g = nt, the equations of motion for the 
minor planet are 


dg 2 d 8 


d 2 /3 

dg 2 


+ 2 


da 

dg 


a, 


p 


3 


In terms of the present coordinates 


n = 


\ (“) + j a2 ‘ f ral ^ ^ ^ 


where 


and 


a 

P 


— - /(I + a) 2 + /3 2 
a 


(1 + a) 2 + /3 2 + - 2 [(1 f cc) cos -£# - p si 


'i -1/2 

n 'iO] j 


p being the distance from the minor planet to Jupiter. Thus 


T2 a - (1 + a) 


-(t ) 3 


- 3a 


+ m 


r a' 
,P*/ 


cos -C# ~ (1 + a.) 


cos tO 


and 


% ~ P 


(if 


+ rn 


_ 

,P> 


sin -td + p 


\—-\ sin 'td 


Later we will also need expressions for the partial derivatives of fi a and fig with 
respect to a, ft, m' and v. These are 


a 


aa 


= - 2 - (|) 3 + 3 (1 + a) 2 (±j + «> /- gj + 3 gj 


cos t0 - (1 + a) 


J 


'} 


^a/3 ~ ^J3 a “ 3 (1 + a) /3 (-pj - 3m' {—) 

\ r j 


a 


cos 'tO - (1 + a) 


sin 'tO + ft 
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n » = 1 -(- 7 ) + 3 / 32 (f) +m ' {"(7) + 3 (f) [(f) sin ' tS + , ' S ] 2 } 


Cl av ~ f "a | 3 (p") ^ + a ) cos ^^“^ sin ^ (-J-) (~) co s ^0 - (1 + a) 


>)’ * > (ffl 


cos 


6 


Cl, 


7 3 V 


= ^3 (1 +a)cos^0-/3sin / t(9 - (~") s * n ^ @ + P 


+ 


(j) * 2 &)] sinie } 


ci , = 1 — n 

am 2 1 + w' av 1 Vp I 


^-^-jcos 't 6 - (1 + a) - cos I 0 


V = 1 n --- - 


2 1 + m' 


> 


(7)T(T) sin ^ e+/? l '(t - ) 2 sin,te 


The dependence on v and part of the dependence on m' comes from the relation 

— = (1 + m') 1/3 v 2/3 . 
a 


The argument, 0 , is treated as an independent quantity, because these expres- 
sions are used for computing the variations of the coefficients in the series 
expansions. 
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For the symmetric periodic orbits, a and Cl a are even functions of 6 While 
fi and are odd so that 


-L 

k=0 


cos k 0 


(X) 


sin k 6 


k=l 


00 


k=0 


fl_= 7 c k cos k 6) 


00 


s k sin k 0 


k=l 


for non-collision orbits* All such series are truncated in the actual computations. 
Putting 


N k =k 


N 


n 


and substituting the series expressions into the equations of motion gives 


(3 + N*) a* + 2N k ^ = - c k 

=-s k J 


\ k = 0, 1, 2, . . 


(A) 


by equating coefficients. If c k and s k were known, the solution would be obtained 
from the formulas 


6 


i 




a 0 


<v 3 




°k = - 


N k 2 -1 


c k + 


Ac = 


N k (N k 2 - 1) 
3 +N 2 


>. k = 1, 2, 3, 


c k - 


N k (N 2 -l) N 2 (N 2 -1) 




(B) 


Since Q a and Hg depend on a and fi , the coefficients c k and s k are not known 
beforehand. However, in many practical cases the solution may be obtained by 
the iterative general perturbations technique. This consists of starting with 
a = /3 = 0, for example, computing approximate values for and s k by har- 
monic analysis of the expressions for Cl a and f <8 , and using these values to solve 
for approximate values of and /3 k . These approximations for cl and /3 are 
then used in a new harmonic analysis of fi a and to give the next approximations 
for and etc. 

This simple process fails in some very interesting cases such as resonance. 
For example, when 


n 


n 



for integer values of p and q , and 



- q 


the formula gives 


N 

p 


= 1 . 


In this case the formulas for Oj, and /3 k have zero divisors at k = p, so the 
method doesn’t work, although the original equations could have a solution with 
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The simple iteration process has been modified to handle cases such ag the 
one just described by taking into account the partial derivatives of c k and s k with 
respect to a and /3. and then computing corrections to an approximate solution. 
This modified process yielded, for example, a linearly stable orbit for n/n' - 
3/2 (the Hilda resonance). This particular orbit will be discussed later. 

The partial derivatives of c k and s k with respect to a. and /3. are obtained 
from the coefficients in the expansions of £l aa , 0. af3 and Q . These functions 
can be represented in the form 


on 


.0 


aa 


- y ' x k cos k 


k=0 


00 


a 


af) 


L 


y,. sin k 0 


k = l 


v* 

-L - 


- ) z,. cos k tJ 


k=0 


Numerical values of the coefficients are computed by harmonic analysis. Then 


3c k 

x k-j * 

da. 

) 

2 


“^,j y| k _ j i +y k+j 

3/3, ' 

2 

3s. 

k _ 

j k - j j +yk + i 

da. 

2 
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where 


3s k 

BA 


lk_j 


k+ j 


('l for k >. j 
l^-l for k < j . 


These variations are the quantities needed for improving the iteration formulas 
which can now be written as 


3 c \ Be 

(3 +N2)S 0lt + 2N k 84 + TJl 8 a, + — 8/3. =e 

L — 1 / — 1 7)/? 


2N k H + N fc 2 84 + 


1=0 1 


r-r* BA 

1=1 'i 


> k=0, 1, 2, . . . (C) 


y 

Z— j Ba. 
i=o 1 


i = 1 


ia. + Y !fl s/3. 
1 ^ ByS. 


'k 


where 


e k = - (3 + N?) a, - 2 N k 4 - c k 


s k = - 2 n k “k - K 4 - s * 


k = 0, 1, 2, 


The quantities So^ and S/ 3 k to be computed are corrections to the approximate 
values 0^ and fi k . 

In the applications the series appearing in this system of equations can be 
truncated at a relatively low order depending on the rate of convergence for the 
particular orbit being computed. The remaining equations are then uncoupled 
in pairs as before. 
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One must also consider the possibility that the symmetric matrix associated 
with this system of linear equations may be singular. The places where this has 
occurred in the applications is at those points on a natural family of orbits where 
the period becomes stationary. These cases occur near a resonance and are 
easily treated by holding the eccentricity fixed rather than the period as in the 
following discussion. 


THE ECCENTRICITY OF THE ORBITS 

In many cases, such as resonance, the eccentricity is a convenient parameter 
to be used in following a family of orbits. For two-body motion the position 
vector may be written as 


r = a P(cos E - e) + a/ 1 e 2 Q sin E 


where P and Q are the usual unit vectors, a is the semi-major axis, e is the 
eccentricity, and E is the eccentric anomaly. The deviations from circular 
motion may be expressed in a and ft by considering 


where 


r - (1 u a) r 0 4 « 


r Q = aP cos g + a Q sin g 
w = -aP sin g + aQcos g 

V 

with g as the mean anomaly. Equating the coefficients of P in the two expressions 
for r , it follows that 


cos E - e = (1 + a) cos g - ft sin g. 

When cos E is expanded in a cosine series in g, the constant term is -e/2. Now 



k=0 


a k cos k 0 
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/? k s i n k 8 , 


k= 1 


and, for the motion to be periodic in the rotating system 

*' « 

6 = g/p 


for some positive integer p . Therefore 


e = (^p " a P )/ 3 * 


This two-body formula may be used to define the mean eccentricity of a perturbed 
orbit. For elliptic motion all the coefficients in the expansions of a and /3 can 
be computed from e using the Bessel functions, and in many cases this is a good 
initial approximation for the perturbed periodic orbit. 

The eccentricity may now be prescribed, and this eliminates one equation in 
the system by putting 


The subscript is placed on e to identify the coefficients with which it is associated. 
The equation which has been eliminated is replaced by another for determining 
the period or the mass ratio as in the following section. 


THE GENERAL PREDICTOR-CORRECTOR FORMULAS 

The Equations (A) relating the coefficients can be written in matrix form as 


DX = f (X) 

where X is the vector of coefficients and A , f (X) is the vector of coefficients 
-c k and -s k which were seen to depend on X , and D is the matrix of coefficients 
multiplying and in Equations (A). It is assumed that the sequence of 
equations is truncated at some appropriate value of k. Note that the matrix 
D depends only on the ratio, v = n/n' , of mean motions for fixed 't . However, 
f depends on m' as well as v . Allowing first order corrections, the equations 
become 
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The components of Bf/Bv and Bf /Bm' are simply the coefficients in the trigonometric 
expansions of Q av , fig v , fi am , and . With v and m' fixed, the iteration 
(or corrector) formulas become 


ax r: - -l^ 1 (f - DX). (D) 

From an established orbit, the predictor formulas which give the variations of 
X with respect to the parameters are 


BX _ / D Bf\ _1 /Bf_ __ BD 
Bv \ BX / \ du Bt- 

and 

BX_ - /p - 1LV 1 — 
Bm 1 \ / Bm' 


(E) 


(F) 


so that the same symmetric matrix is to be inverted in each case. 

The eccentricity, e p , is introduced as a parameter using the relation 


S/3 = 3 S e + Sct 

p p 

wherever S/3 p occurs in the equations. Then, when Se p is specified, either Sv . 
or Sm' is obtained as part of the solution. This is accomplished by simple 
manipulations of the vectors and matrices in the above formulas. After the 
solution, the new value of /3 p is computed from 


A, =3 


+ a 


12 


mm* 


mmm 
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THE LINEAR STABILITY ANALYSIS 


Putting 


7 = da/dg 


and 


8 = dyd/dg 

the equations of motion can be written as 


da/dg = 7 
d/S/dg = S 

d 7 / dg = 3 a + 2 8 + Cl a 
dS/dg = - 27 + Hq 

so that the variational equations are 


with 


d<&/ dg = A<E> 


0 0 10 


A = 


0 

3 + fi 

CLCL 


0 0 1 

% 0 2 




- 2 



The functions H , fi 0 and 0. oo are available in series form from the 

CLCL CLp p p 

iteration process, so the variational equations can be integrated directly. For 
the present this is being done numerically using power series. 


THE APPLICATIONS 

Previous applications were mentioned in the introduction. Some additional 
orbits and families of orbits are given here as illustrations. These orbits all 
« close after one revolution around the sun in the synodic system (•£ = 1 in the 
definition of the trigonometric argument). 

/ 
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A typical stable orbit of the first kind is shown in the synodic coordinate 
system in Figure la. Such orbits are nearly circular and very easy to compute. 
The deviations from the circular reference motion are shown in the a , ft plot of 
Figure lb, and the series coefficients are given in Table 1. The values are given 
to six decimals for purposes of illustration. The orbits were all computed to an 
accuracy of about 12 decimals. 

As the ratio of mean motions, v - n /n' , is decreased toward the resonance 
value, v = 2, the deviations from circular motion increase as shown in Figure 2a 
for the stable orbit at v = 2.001. The a, ft plot of Figure 2b for this orbit con- 
sists of two loops which are nearly identical. The series coefficients, given in 
Table 2, are much larger now, but the rate of convergence is still good. This 
orbit is best described as a slowly rotating perturbed ellipse with eccentricity 
e 2 = 0.388739. The coefficients with odd subscripts would all be zero if the orbit 
were a true ellipse. 

There are more periodic orbits which are nearly circular for 3/2 < v < 2/1. 
The deviations from circular motion for one such orbit are shown in Figure 3 
where v = 1.6. This orbit is also stable and has the series coefficients shown 
in Table 3. 

As the ratio of mean motions is reduced toward the 3/2 resonance value, the 
eccentricity increases again. There is a stable orbit at the exact resonance 
shown in Figure 4a. Now the a, ft plot, Figure 4b, consists of three nearly 
identical loops. The dominant coefficients (with subscripts which are multiples 
of 3) are given in Table 4. This orbit is a slightly perturbed ellipse with e 3 = 
0.453692 and no secular motion of the perihelion. 

Each of the above orbits was computed as a member of a natural family 
with m' = 1/1047.35 and using v = n /n' as the parameter. For the orbit of 
Figure 4 this would not have been possible without using the partial derivatives 
in the iteration because of the zero divisors. For the orbits near the resonances, 
v = 2/1 and v =3/2, the eccentricity, e 2 or e 3 , would serve as a better 
parameter. 

Each of the orbits except the one of Figure 4 has also been computed starting 
from circular two-body motion (m' =0) and then increasing m' to the value for 
Jupiter while holding v fixed. This is the technique for orbits of the first kind 
(Poincare). The variations of e 2 with m' for three families of this type are 
shown in Figure 5. 

Each of the orbits has also been computed by starting from elliptic two-body 
motion (m' =0) and then increasing m' to the value for Jupiter while holding the 
eccentricity e p fixed as in the method for orbits of the second kind (Poincare"). 
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In this approach p = 2 for the orbits of Figures 1 and 2 while p = 3 for the orbit 
of Figure 4. As a numerical experiment the orbit of Figure 3 was computed by 
continuation from two different two-body elliptic orbits, one starting from v - 2 
and holding e 2 fixed, and the other starting from v = 3/2 and holding e 3 fixed. 

The variations of v with respect to m' for these two families along with three 
others are shown in Figure 6. 

In following families of the types illustrated in Figures 5 and 6 there will be 
cases where the curves have vertical tangents (the mass m' is stationary at a 
point on the family). The mass m' is not a suitable parameter in a neighborhood 
of such a point, and the matrix involved in the computation becomes singular. 

In such cases the ratio of mean motions, v , and the eccentricity, e p , are 
specified, while m' is obtained in the solution along with the other coefficients. 

This is somewhat akin to prescribing a function and then searching for the 
problem which it solves, but in this case the procedure is necessary and justified 
for the purpose of continuing a family of orbits. 

The relationship between the a , {3 and synodic coordinate systems is indicated 
in Figure 7. 


15 








y 



Figure la-Nearly circular stable periodic orbit of the finst hand with mean motion slightly greater 
than that of Hecuba resonance, m ' = 1/1047.35, v-n/n' ~ 2.1, - 0.011905. (Synodic coor- 

dinate system with the sun at the origin and Jupiter at X = 1.) 
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Figure lb-Deviations from circular motion for the orbit of Figure la. The two loops cross at 
the top of the figure, while at the bottom they are nearly tangent but do not cross. Uniform 
circular motion, would give a = /3 = 0 throughout the orbit. The sun is at a = - 1 while 
Jupiter moves on the circle (1 + a.) 2 + /3 2 = (a 1 /a) 2 . 


Table 1 


Coefficients of the trigonometric series 
for a and j3 in the orbit of Figures la and 
lb 


k 

\ *io 6 

4-io 6 

0 

-122 

0 

1 

595 

-2028 

2 

-12158 

23556 

3 

-491 

719 

4 

-60 

202 

5 

-43 

57 

6 

-18 

23 

7 

-8 

10 

8 

-4 

4 

9 

-2 

2 


-1 

1 

11 

-1 

1 
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Table 2 

Coefficients of the trigonometric series 
for the orbit of Figures 2a and 2b 


k <vl0 6 /\.10 6 


0 -75926 0 

1 -1433 2179 

2 -410408 755808 

3 +198 334 

4 68305 28584 

5 +149 128 

6 18652 13566 

7 +72 65 

8 6429 5342 

9 35 33 

10 2429 2137 

11 17 16 

12 968 879 

13 8 8 

14 400 370 

15 4 4 

16 170 159 

17 2 2 

18 73 69 

19 1 1 

20 32 31 

22 14 14 

24 6 6 

26 3 3 

28 1 1 

30 1 1 








Figure 3— Deviations from near circular motion of a stable periodic orbit between the 2/1 and 3/2 
resonances, m' = 1/1047.35, v = 1 .6, e 2 =-0.008 292, e 3 = 0.007 791. 
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Table 3 

Coefficients of the trigonometric series 
for the orbit of Figure 3 


k 

a, -10 6 

k 

A,-io 6 

0 

-235 

0 

1 

1842 

-8567 

2 

7202 

-17674 

3 

-8168 

15207 

4 

-1234 

1935 

5 

-522 

585 

6 

-185 

271 

7 

-100 

128 

8 

-56 

66 

9 

-32 

37 

10 

-19 

21 

11 

-11 

13 

12 

-7 

8 

13 

-4 

5 

14 

-3 

3 

15 

-2 

2 

16 

-1 

1 

17 

-1 

1 

18 

-1 

1 
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Figure 4b-Deviations from circular motion for the orbit of Figure 4a. The figure consists of 
three nearly identical loops similar to those of Figure 2b except that the eccentricity is now 






Table 4 

Dominant coefficients of the trigonometric 
series for the orbit of Figures 4a and 4b 


k 

<vio 6 

£ k -io 6 

0 

-103710 

0 

3 

-487824 

873251 

6 

89846 

32529 

9 

27907 

19671 

12 

11018 

8978 

15 

4778 

4144 

18 

2187 

1963 

21 

1038 

952 

24 

506 

471 

27 

252 

236 

30 

127 

120 

33 

65 

62 

36 

34 

32 

39 

18 

17 

42 

9 

9 

45 

5 

5 

48 

3 

3 

51 

1 

1 

54 

1 

1 
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z/=2.001 


z/=2.01 


^= 1.9835 


0.0002 


0.0004 


0.0006 


0.0008 


Figure 5-Three families of the first kind. These families start from circular motion at m ' - 0 
and are continued up to m* = 1/1047.35 holding the ratio of mean motions, V- n/n*, at a 
fixed value. These curves show the variations of the eccentricity ej for three values of v 
near 2. The cases v = p/(p— 1) for integral p > 1 are singular for continuations of this type. 
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e 2 =0.0005 


e 2 =o.ooi 



e 2 =-0. 008292 


e 3 =+0.007791 


0.0002 


0.0004 


0.0006 


0.0008 


Figure 6-Five families of the second kind. These families start from elliptic motion at exact 
resonance v=p/(p-l) form' =0 and are continued up to m 1 = 1 /1047. 35 holding the eccen- 
tricity e p fixed. The ratio of mean motions, v, is computed as part of the solution. Curves are 
shown for four families starting at v - 2 and one family starting at v = 3/2. The families e ^ 
-0.008 292 and e^ = 0.007 791 intersect at the orbit of Figure 3. The cases e p = 0 are singu- 
lar for continuations of this type from v- p/(p-l). 




SUN 


r 


JUPITER 


Figure 7— Geometry of the coordinate systems. 

|r Q | = a, |w| = a, |r| = a [(1 + a) 2 + /3 2 ] 1,/2 > |r # l - a ‘ , % & - (n-n')t, 

I P I = a i"( 1 + a ) 2 + ft 2 + - 2 [(1 + a ) cos^6 , -/3sinT<9j'l 




CONTROLS ON THE COMPUTATION 


The Jacobi constant plays no role in the computational scheme used in this 
study. Its value is monitored as an indicator of numerical inaccuracies or 
difficulties. 

The iteration process is continued until the residuals in the differential 
equations fall below an acceptable tolerance. The number of terms kept in the 
series and the order of the matrix used in the iteration are determined for each 
case so as to give the required accuracy and convergence. Generally it was 
required that the last ten terms computed in the series be less than 0.5*10 ' 12 . 

When the iteration matrix becomes ill-conditioned, a change is made to a 
different independent parameter. No cases have been encountered where this 
presents jmy serious difficulty. 

Several cases have been checked by comparison with numerical integration, 
and the expressions for the variations have been checked by numerical differencing 
of the series coefficients. 


COMPARISON WITH NUMERICAL INTEGRATION 

The author is aware of no numerical integration program which has the 
facility for continuing families of orbits in all of the ways discussed here, or 
which gives the variations with respect to a parameter of a representation of 
the entire orbit such as those given in equations (E) and (F). On the other hand, 
the present technique is not suitable when the series convergence becomes too 
slow, nor is it useful for computing orbits which are not periodic. The necessary 
modifications for regularization have not yet been made. 

It should be mentioned that, although the restricted problem of three bodies 
is extremely rich in natural families (Wintner) of periodic orbits, elements of 
many of these families have been discovered only by accident. Continuation from 
two-body motion as suggested by Poincare" and applied here provides a basis for 
the systematic exploration of many of these families of orbits. 

The present method can be used in conjunction with a good numerical inte- 
gration program taking full advantage of the special features of each. 
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CONCLUSION 


The present study demonstrates the strength of a semi-analytic development 
in cases where small parameter methods fail. 

The method described here has been justified mathematically for isolated 
periodic solutions (see Urabe (1965) and Stokes (1969)), but further work is 
needed for singular solutions (periodic orbits which are members of a natural 
family) as in this study. The indeterminancies are removed here by applying 
constraints. The value of a suitable parameter is specified, and the epoch of 
time is defined to occur at a perpendicular crossing of the axis for symmetric 
orbits. For non-symmetric orbits one might require one of the leading coefficients 
in the series to be zero to tie down the epoch. 

In the restricted problem of four-bodies, once the motions of the three pri- 
maries are determined, the periodic solutions for the infinitesimal body are 
isolated, so these difficulties are not encountered. Nor would they be present in 
the reduced (or elliptic) problem of three-bodies. 

There are many possible applications of this method in celestial mechanics 
problems, and several are being explored. 
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